This files shows the basic steps of fly eye measurement and feature extraction with one image as an example.
# Load packages
library(lattice)
library(ggplot2)
library(sp)
library(raster)
library(tiff)
library(EBImage)
library(Gmedian)
library(gridExtra)
# source helpful R code and specify image folder path
source('./src_Fly.eye.pat_funcs.R')
dir.in='./example_images/'
dir.out='./output/'
cat("check dir: ",dir.in,"; dirout",dir.out,"\n");
## check dir: ./example_images/ ; dirout ./output/
# Read in images and begin process
image.filenames <- list.files(path=dir.in,pattern="*jpg$", full.name=F)
image.filenames=image.filenames[1]; #use one image as an example
print(image.filenames);
## [1] "test1.jpg"
n.images=length(image.filenames);
image.files=paste(dir.in,image.filenames,sep='/');
# process one image at a time
# save each image quantitive features in a lsit
image.i = 1
image = readImage(image.files[[image.i]] )
image.name = image.filenames[image.i]
cat('begin processing image',image.name,',',image.i,'images out of',n.images,'total images\n')
## begin processing image test1.jpg , 1 images out of 1 total images
display(image, method="raster")
dim(image)
## [1] 2560 1920 3
# Resize to fit memory
image.res=resFunc(image)
display(image.res, method="raster")
dim(image.res)
## [1] 640 480 3
# Assign resized images RGB channels to data frames
image.Ori <- RGBintoDF(image.res)
head(image.Ori)
## x y R G B
## 1 1 480 1 1 1
## 2 2 480 1 1 1
## 3 3 480 1 1 1
## 4 4 480 1 1 1
## 5 5 480 1 1 1
## 6 6 480 1 1 1
# White TopHat morphological transform
image.Top <- wTopHat(image.res,y=5,z='diamond')
# Display images
par(mfcol=c(1,2))
dispImg(image.res)
dispImgT(image.Top,0.99)
## Threshold and centroids
# Assign gray channel to data frame
image.grey <- GintoDF(image.Top)
head(image.grey)
## x y G
## 1 1 480 0
## 2 2 480 0
## 3 3 480 0
## 4 4 480 0
## 5 5 480 0
## 6 6 480 0
# Threshold to keep pixels with intensity > 0.99 quantile
image.thres.retain<-image.grey[image.grey$G > quantile(image.grey$G,0.99), ]
# select ROI, region of interest through two rounds of pixel fitering based on distances to centroids
image.thres = image.thres.retain
cutoffs=c(0.8,0.3);
for(cutoff in cutoffs){
cat('ROI selection with cutoff',cutoff,'\n')
# Estimate images centroid
centroids <- Weiszfeld(image.thres[,1:2])
# extract each pixel x,y coordiantes retained pixels
thresXY <- image.thres[,1:2,drop=FALSE]
p1<-ggplot(thresXY,aes(x=x,y=y))+geom_point()+
geom_point(aes(x=centroids$median[1], y=centroids$median[2]),colour="red",size=4)+
theme_bw()+ ggtitle(image.name)+xlab('')+ylab('')+
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank())
## Distances to centroide,calculate distances to centroid
pdist <- pointDistance(p1=thresXY, p2=centroids$median,lonlat=F)
# Mark distances > 0.8 quantile, as they belong to
# points outside the eye boundary in their majority
pLogic <- (pdist < quantile(pdist, cutoff));
distCent <- cbind(distCent = pdist, selected = pLogic)
# Plot example histograms
p2<-ggplot(data.frame(distCent),aes(x=distCent))+geom_histogram(bins = 30)+
geom_vline(xintercept = quantile(distCent[,1], cutoff),col="red",lty="dashed", lwd=2)+
theme_bw()+ ggtitle(image.name)
# Join thresholded and distances lists
thresDist <- cbind(image.thres,distCent);
# retain points with distance < quantile cutoff 0.8
distSelect <- thresDist[thresDist$selected==1,]
# Plot examples (black clouds with blue roi overlay)
p3 <- ggplot(data=thresDist, aes(x=x, y=y,color=selected))+
geom_point(show.legend = FALSE) +
theme_bw()+ ggtitle(image.name)+xlab('')+ylab('')+
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank())
print(grid.arrange(p1,p2,p3,ncol=3))
# Overlay to image and plot example
p <- ggplot(data = image.Ori, aes(x = x, y = y)) +
geom_point(colour = rgb(image.Ori[c("R","G","B")])) +
labs(title = "Original Eye selected Points",
cex=0.5) +xlab("x") +ylab("y") +
geom_point(data=distSelect, alpha=0.2) +
plotTheme()
image.thres <- distSelect[,1:3]
}
## ROI selection with cutoff 0.8
## TableGrob (1 x 3) "arrange": 3 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (1-1,3-3) arrange gtable[layout]
## ROI selection with cutoff 0.3
## TableGrob (1 x 3) "arrange": 3 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (1-1,3-3) arrange gtable[layout]
## Subset and confidence ellipse and Add ellipse to plot
ellPlots = ellPlot(distSelect,0.90)
ellPlots #it's a ggplot object
# Extract components
build = ggplot_build(ellPlots)$data
ell = build[[2]]
# Select original image points inside ellipse
orig.ell = data.frame(image.Ori[,1:5],
in.ell = as.logical(point.in.polygon(image.Ori$x,image.Ori$y, ell$x, ell$y)))
orig.pix = orig.ell[orig.ell$in.ell==TRUE,]
# plot example: raw eye image + ROI circled in blue
p <- ggplot(data = image.Ori, aes(x = x, y = y)) +
geom_point(colour = rgb(image.Ori[c("R", "G","B")])) +
labs(title = "Original Eye final ROI", cex=0.5) +
xlab("x") + ylab("y") +
geom_polygon(data=ell[,1:2], alpha=0.2,
size=1, color="blue") +#plotTheme()
theme_bw()+xlab('')+ylab('')+
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank())
print(p)
## Create image from final ROI
rois.image = roitoImg(orig.pix)
display(rois.image,method='raster'); #check image
## begin extract ommatidia configuration features
pic=rois.image
#hist(pic);grid()
y = equalize(pic)
#hist(y);grid()
#display(y, title='Equalized Grayscale Image',method='raster')
grayimage<-channel(y,"grey")
#display(grayimage,method='raster')
nmask = thresh(grayimage, w=5, h=5, offset=0.02); #display(nmask)
nmask = opening(nmask, makeBrush(3, shape='disc')); #display(nmask)
nmask = fillHull(nmask); #display(nmask)
nmask = bwlabel(nmask)
cat("Number of ommatidia,",max(nmask),"\n");
## Number of ommatidia, 75
fts = data.frame(computeFeatures.moment(nmask))
head(fts)
## m.cx m.cy m.majoraxis m.eccentricity m.theta
## 1 50.50000 1.500000 4.472136 0.8944272 0.0000000
## 2 38.06452 5.903226 8.830081 0.8437787 1.5130830
## 3 62.42623 5.459016 18.670999 0.9729211 0.3624306
## 4 45.89474 6.421053 5.544698 0.6322809 -0.2237600
## 5 31.25926 9.925926 6.518058 0.5827875 -0.1450000
## 6 54.40741 10.925926 6.421501 0.5724126 0.2180192
fts2 <- data.frame(computeFeatures.shape(nmask))
head(fts2)
## s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
## 1 8 8 1.144123 0.4370160 0.7071068 1.581139
## 2 31 19 2.759500 0.8339106 1.1263650 3.846429
## 3 61 35 4.752899 2.2218203 1.5452363 8.321352
## 4 19 13 1.986977 0.3333103 1.3931362 2.560510
## 5 27 16 2.465936 0.3606885 1.9465514 3.006503
## 6 27 15 2.471452 0.3104342 1.8961950 2.988125
par(mfrow=c(1,5));
display(pic,"raster")
display(y, title='Equalized Grayscale Image',method='raster')
display(grayimage,method='raster')
display(nmask,"raster");
display(nmask,"raster");
text(fts[,"m.cx"], fts[,"m.cy"],
labels=seq_len(nrow(fts)), col="red", cex=0.8)
par(mfrow=c(1,2));
display(pic,"raster")
#display(y, title='Equalized Grayscale Image',method='raster')
#display(grayimage,method='raster')
display(nmask,"raster");
#display(nmask,"raster");
#text(fts[,"m.cx"], fts[,"m.cy"],labels=seq_len(nrow(fts)), col="red", cex=0.8)
# fts and fts2 contain raw measurements on each ommatidium, which can be written out to a file
#out.file1=paste(dir.out,"/",image.name,"-coord.txt",sep='');
#out.file2=paste(dir.out,"/",image.name,"-area.txt",sep='');
#write.table(fts,file=out.file1,quote=F,row.names = F);
#write.table(fts2,file=out.file2,quote=F,row.names = F);
## ROI seleciton done, begin to compute features based on raw ommadium measurements
## some filter:
# images contain less then 6 ommaditia, discard
# when calculating each feature, remove some extreme values, the largest 3 and smallest 3.
# number of ommatida
n.ommatidia=nrow(fts)
if(n.ommatidia<6){ result[[image.name]]=NA;next}
# There are four features computed from data.frame fts
# nn: pairwise ommatidia nearest neighbor distances
# cc: each ommatidium essentricity
# nn.mean, nn.sd, cc.mean, cc.sd
colnames(fts)[c(1,2,4)]; #x,y coordinates, essentricity
## [1] "m.cx" "m.cy" "m.eccentricity"
nn.dist=c();
for(i in 1:nrow(fts)){
pairwise.dist=c();
for(j in 1:nrow(fts)){
if(i==j){next}
distance=((fts[i,]$m.cx-fts[j,]$m.cx)^2+(fts[i,]$m.cy-fts[j,]$m.cy)^2)^0.05;
pairwise.dist=c(pairwise.dist,distance);
}
nn.dist=c(nn.dist,min(pairwise.dist))
}
#length(nn.dist)==nrow(fts)
nn.out=get_mean_sd(nn.dist)
names(nn.out)=c('nn.mean','nn.sd')
cc.out=get_mean_sd(fts$m.eccentricity)
names(cc.out)=c('cc.mean','cc.sd');
x0=c(n.ommatidia,nn.out,cc.out)
names(x0)=c('n.ommatidia','nn.mean','nn.sd','cc.mean','cc.sd')
# There are 12 features computed from data.frame fts2
# mean and sd for each column variable
colnames(fts2)
## [1] "s.area" "s.perimeter" "s.radius.mean" "s.radius.sd"
## [5] "s.radius.min" "s.radius.max"
out=apply(fts2,2,get_mean_sd)
x1=out[1,]
x2=out[2,]
name1=paste0('mean(',names(x1),')')
name2=paste0('sd(',names(x1),')')
names(x1)=name1;names(x2)=name2
out=c(x0,x1,x2)
plot.new()
grid.table(data.frame(out))
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] gridExtra_2.3 Gmedian_1.2.4 EBImage_4.28.0 tiff_0.1-5
## [5] raster_3.0-7 sp_1.3-2 ggplot2_3.3.2 lattice_0.20-38
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.0 locfit_1.5-9.1 xfun_0.11
## [4] purrr_0.3.3 colorspace_1.4-1 vctrs_0.3.2
## [7] generics_0.0.2 htmltools_0.4.0 yaml_2.2.0
## [10] rlang_0.4.7 pillar_1.4.6 glue_1.4.2
## [13] withr_2.1.2 BiocGenerics_0.32.0 jpeg_0.1-8.1
## [16] lifecycle_0.2.0 robustbase_0.93-5 stringr_1.4.0
## [19] munsell_0.5.0 gtable_0.3.0 htmlwidgets_1.5.1
## [22] codetools_0.2-16 evaluate_0.14 labeling_0.3
## [25] knitr_1.26 parallel_3.6.1 DEoptimR_1.0-8
## [28] Rcpp_1.0.3 scales_1.1.1 abind_1.4-5
## [31] farver_2.0.3 RSpectra_0.15-0 png_0.1-7
## [34] digest_0.6.22 stringi_1.4.3 dplyr_1.0.2
## [37] grid_3.6.1 tools_3.6.1 bitops_1.0-6
## [40] magrittr_1.5 RCurl_1.95-4.12 tibble_3.0.3
## [43] crayon_1.3.4 pkgconfig_2.0.3 ellipsis_0.3.0
## [46] MASS_7.3-51.4 Matrix_1.2-17 rmarkdown_1.17
## [49] R6_2.4.1 fftwtools_0.9-8 compiler_3.6.1
#installed.packages()[names(sessionInfo()$otherPkgs), "Version"]